## Warning: package 'dplyr' was built under R version 3.5.2
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

1 Goal

This notebook can be called with knit_notebook8_byESM.R; it can be called in a loop over ESMs and target experiments, and will knit an analysis html for each.

These values must be defined globally in knit_notebook8_byESM.R or uncommented and defined here to run this notebook directly:

# # what we will emulate:
# esm_name <- "MPI-ESM1-2-LR"
# esm_experiment <- 'ssp245'
# 
# # How many times do we want to draw full generated ensembles:
# Ndraws <- 1

print(esm_name)
## [1] "CanESM5"
print(esm_experiment)
## [1] "ssp245"

We will only explore tol=0.1 degC for matching in this notebook (per CT, from the SD of the smoothed Tgav time series. Should check for each ESM, for now just use 0.1 value from CanESM).

2 Set Up

# Then load the functions we will want to use, these are currently written in R and will be translated into python. 
source("functions/nearest_neighbor_matching.R") # the function match_neighborhood is defined here
source("functions/permuting_matches.R")  # the function permute_stitching_recipe is defined here
source("functions/stitching_functions.R")  # the function stitch_global_mean is definend here 
source("functions/gridded_recipe.R") # wrappers for the work flow from matching to output csv needed by python for 
                                     # layering in gridded data to output as netcdf

set.seed(1)

Load the inputs that will be used. In addition to SSP534-over having enough data issues to fully exclude, at least some of the realizations for SSP434 and SSP460 just don’t have future data at all for CanESM5. And for other models, some realizations have odd windows (Like a 2012-2014 window) and those realizations are excluded as well. We need to go back to pangeo and fix this for sure, but for now, I’ll work with the archive I have. Because we don’t want to lose too much of the archive now, I will only remove those realizations and not a full experiment for a single bad data point.

# Time series of the raw global mean temperature anomaly, this is the data that is 
# going to be stitched together. Historical data is pasted into every scenario. No
# smoothing has been done.
tgav_data <- read.csv("inputs/main_raw_pasted_tgav_anomaly_all_pangeo_list_models.csv", 
                      stringsAsFactors = FALSE) %>% dplyr::select(-X)


# A chunked smoothed tgav anomaly archive of data. This was tgav_data but several 
# steps were taken in python to transform it into the "chunked" data and save it 
# so that we do not have to repeate this process so many times. 
read.csv(paste0('inputs/', esm_name, '_archive_data.csv'), stringsAsFactors = FALSE) %>%
  # Exclude ssp534-over:
  dplyr::filter(experiment != 'ssp534-over') %>%
  # eliminate any individual realizationXexperiment runs that don't have all the years:
  group_by(experiment, variable,   model  , ensemble) %>%
  mutate(nWindows = n()) %>%
  ungroup %>%
  filter(nWindows == 28) %>%
  select(-nWindows) ->
  archive_data


# The target data - all ssp245 realizations (smoothed because it's the target for
# matching)
archive_data %>%
  filter(experiment == esm_experiment) ->
  target_data


# The corresponding unsmoothed original data for validation
tgav_data %>%  
  filter(model == unique(target_data$model) & experiment == unique(target_data$experiment))  %>%
  select(-file, -timestep, -grid_type) -> 
  original_data

We also need to prepare the file of path metadata for pangeo. This is a csv file that contains the pangeo paths to different variables of interest. By selecting the paths for tas, psl, and pr, keeping only those experimentXensemble combinations that have defined paths to netcdfs for all three variables, and then using the resulting experiment X ensemble combinations as an additional filter on the archives we use for matching, we make sure that we can then stitch together gridded data for all 3 of these variables no matter what our generated ensembles are. Adding or removing variables of interest is relatively trivial. For now, the path metadata file has information for tas, psl, pr, tasmin, and tasmax.

# The main pangeo file of path metadata so that we can filter the archives to include only
# the experiment*ensemble members with netcdfs for all variables of interest
read.csv("inputs/pangeo_path_metadata_tas_psl_pr_tmax_tmin.csv", 
         stringsAsFactors = FALSE,
         na.strings = c("", " ", "NA")) %>%  
  rename(model = source_id, experiment = experiment_id, ensemble = member_id) %>%
  filter(model == unique(target_data$model) )  %>%
  select(model, experiment, ensemble, table_id, tas_zstore, psl_zstore, pr_zstore) %>%
  # cut out any rows that don't have paths to netcdfs for all variables:
  na.omit %>%
  select(model, experiment, ensemble)  %>%
  unite("acceptable", c(model, experiment, ensemble), sep="~")->
  acceptable_archive_entries_from_path_metadata

3 Generate the new ensembles

We’re actually going to generate multiple different ensembles for a target of ssp245, based on two different archives.

Archive 1 will exclude the ssp245 runs from potential matches.

Archive 2 will include the ssp245 runs - when we get to the point of targeting SSP585 runs, I think we will have to to include those runs in the archive to get at all decent generated Tgavs at the end of the century when temperature is higher. It’s still not just regurgitating the target ensembles, it will include mixing and matching from those parts.

We will also do some light exploration of the matching neighborhood size (tol argument).

# Archive data not including the target experiment:
archive_data %>% 
  # Only match to the same ESM:
  dplyr::filter(model == unique(target_data$model)) %>%  
  # Don't match to any ensemble members in the target experiment:
  dplyr::filter(experiment != unique(target_data$experiment) ) %>%
  # Final filter to make sure we only use the experiment*ensemble members with pangeo
  # path metadata for netcdfs for all varaibles of interest
  unite("acceptable", c(model, experiment, ensemble), sep="~", remove=FALSE) %>%
  filter(acceptable %in% unique(acceptable_archive_entries_from_path_metadata$acceptable)) %>%
  select(-acceptable) ->
  archive_subset1 

# what does this archive look like:
archive_subset1 %>% 
  select(experiment, ensemble) %>%
  distinct() %>% 
  group_by(experiment) %>% 
  summarize(n_complete_ensemble_members = n()) %>%
  ungroup %>%
  kable()
experiment n_complete_ensemble_members
ssp119 20
ssp126 25
ssp370 25
ssp434 5
ssp460 5
ssp585 25
# Archive data including the target experiment:
archive_data %>% 
  # Only match to the same ESM:
  dplyr::filter(model == unique(target_data$model)) %>%  
  # Final filter to make sure we only use the experiment*ensemble members with pangeo
  # path metadata for netcdfs for all varaibles of interest
  unite("acceptable", c(model, experiment, ensemble), sep="~", remove=FALSE) %>%
  filter(acceptable %in% unique(acceptable_archive_entries_from_path_metadata$acceptable)) %>%
  select(-acceptable) ->
  archive_subset2 

# what does this archive look like:
archive_subset2 %>% 
  select(experiment, ensemble) %>%
  distinct() %>% 
  group_by(experiment) %>% 
  summarize(n_complete_ensemble_members = n()) %>%
  ungroup %>%
  kable()
experiment n_complete_ensemble_members
ssp119 20
ssp126 25
ssp245 25
ssp370 25
ssp434 5
ssp460 5
ssp585 25
for (draw in 1:Ndraws){
  
  # The matching to the archive without ssp245 results
  match_neighborhood(target_data = target_data,
                     archive_data = archive_subset1,
                     tol=0.1) -> matched_data
  matched_data %>%
    # Convert them to a sample of individual Tgav Recipes:
    permute_stitching_recipes(N_matches = 2000, matched_data = .,
                              archive = archive_subset1) ->
    recipes_1_0p1 
  
  # convert these recipes to something that can be ingested by python 
  # to layer in gridded netcdfs and save
  gridded_recipes_1_0p1 <- generate_gridded_recipe(recipes_1_0p1)
  write.csv(gridded_recipes_1_0p1,
            paste0('generated_ensembles/', esm_name, '/gridded_recipes_for_python_', esm_experiment,'_draw', draw, '_experiment_not_in_matching.csv'),
            row.names = FALSE)
  
  # Proceed to stitch the actual Tgavs for analysis here.
  recipes_1_0p1 %>%
    # And stitch the recipes into tgavs:
    stitch_global_mean(recipes=., data = tgav_data) %>%
    # add some identifying info
    mutate(tol = 0.1,
           archive = paste0('without_', esm_experiment)) ->
    rslts_1_0p1

  # Check for collpase in these results
    rslts_1_0p1 %>%
    select(year, value, stitching_id) %>%
    spread(year, value) ->
    wide

    # Collapse would have occured if any of these trajectories go through the same value
    # in the same year -> if the length of distinct values in any column is less than
    # the number of stitching_id's
    counts <- apply(wide, 2, length)

    print(paste('Archive without', esm_experiment,'; tol = 0.1: Is there any collapse in the generated ensemble?', any(counts != length(unique(wide$stitching_id)))))

    

  # The matching to the archive with target experiment included results
  match_neighborhood(target_data = target_data,
                     archive_data = archive_subset2,
                     tol=0.1) %>%
    # Convert them to a sample of individual Tgav Recipes:
    permute_stitching_recipes(N_matches = 2000, matched_data = .,
                              archive = archive_subset2) ->
    recipes_2_0p1
  
  # convert these recipes to something that can be ingested by python 
  # to layer in gridded netcdfs and save
  gridded_recipes_2_0p1 <- generate_gridded_recipe(recipes_2_0p1)
  write.csv(gridded_recipes_2_0p1,
            paste0('generated_ensembles/', esm_name, '/gridded_recipes_for_python_', esm_experiment,'_draw', draw, '_experiment_in_matching.csv'),
            row.names = FALSE)
  
  
  recipes_2_0p1 %>%
    # And stitch the recipes into tgavs:
    stitch_global_mean(recipes=., data = tgav_data) %>%
    # add some identifying info
    mutate(tol = 0.1,
           archive = paste0('with_', esm_experiment)) ->
    rslts_2_0p1

    # Check for collpase in these results
    rslts_2_0p1 %>%
    select(year, value, stitching_id) %>%
    spread(year, value) ->
    wide

    # Collapse would have occured if any of these trajectories go through the same value
    # in the same year -> if the length of distinct values in any column is less than
    # the number of stitching_id's
    counts <- apply(wide, 2, length)

    print(paste('Archive with', esm_experiment,'; tol = 0.1: Is there any collapse in the generated ensemble?', any(counts != length(unique(wide$stitching_id)))))

    
    
    
  # save it all off  
  write.csv(bind_rows(rslts_1_0p1, 
                      rslts_2_0p1) %>%
              mutate(draw = draw),
            paste0('generated_ensembles/', esm_name, '/generated_ensembles_', esm_experiment,'_draw', draw, '.csv'),
            row.names = FALSE)
}
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp245 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"

4 Validate - ensembles

Comparing the generated ensembles to the training ensemble with error metrics along the lines of Tebaldi et al 2020.

That is, we need the means and variances of both the target ensemble and the generated ensemble. Specifically, mean across time and ensemble members, variance across time and ensemble members.

# values of the target data
mean_target <- mean(original_data$value)
mean_target_hist <- mean(filter(original_data, year <= 2014)$value)
mean_target_future <- mean(filter(original_data, year > 2014)$value)
  
var_target <- sd(original_data$value)
var_target_hist <- sd(filter(original_data, year <= 2014)$value)
var_target_future <- sd(filter(original_data, year > 2014)$value)
  

# generated ensemble files
# read in all generated ensembles.
files <- list.files(path = paste0('generated_ensembles/', esm_name), 
                    pattern = paste0('generated_ensembles_', esm_experiment), 
                    full.names = TRUE) 

for (i in 1:length(files)){

    # Metrics for this draw #################################################################
  generated_ensemble_draw <- read.csv(files[i], stringsAsFactors = FALSE)
  
  # Values of the generated data and then the metrics
  generated_ensemble_draw %>% 
    group_by(tol, archive, draw) %>%
    summarize(mean_gen = mean(value), 
              var_gen = sd(value)) %>%
    ungroup %>%
    mutate(time_period = 'hist+future',
           mean_target = mean_target,
           var_target = var_target,
           bias = abs(mean_target - mean_gen),
           E1 = bias/var_target,
           E2 = var_gen/var_target) ->
    metrics
  
  
  generated_ensemble_draw %>% 
    filter(year <= 2014) %>%
    group_by(tol, archive, draw) %>%
    summarize(mean_gen = mean(value), 
              var_gen = sd(value)) %>%
    ungroup %>%
    mutate(time_period = 'hist',
           mean_target = mean_target_hist,
           var_target = var_target_hist,
           bias = abs(mean_target - mean_gen),
           E1 = bias/var_target,
           E2 = var_gen/var_target) ->
    metrics_hist
  
  
  generated_ensemble_draw%>% 
    filter(year > 2014) %>%
    group_by(tol, archive, draw) %>%
    summarize(mean_gen = mean(value), 
              var_gen = sd(value)) %>%
    ungroup %>%
    mutate(time_period = 'future',
           mean_target = mean_target_future,
           var_target = var_target_future,
           bias = abs(mean_target - mean_gen),
           E1 = bias/var_target,
           E2 = var_gen/var_target) ->
    metrics_future
  
  
  generated_ensemble_draw %>%
    select(stitching_id, tol, archive, draw) %>%
    distinct %>%
    group_by(tol, archive, draw) %>%
    summarise(n_stitched = n()) %>%
    ungroup ->
    tmp
  
  bind_rows(tmp %>%
              left_join(metrics, by = c('tol', 'archive', 'draw')),
            tmp %>%
              left_join(metrics_hist, by = c('tol', 'archive', 'draw')) ,
            tmp %>%
              left_join(metrics_future, by = c('tol', 'archive', 'draw')))  ->
    plotting_tbl
  rm(tmp)
  
   write.csv(plotting_tbl %>%
              mutate(draw = draw),
           paste0('generated_ensembles/', esm_name ,'/metrics_', esm_experiment,'_draw', i, '.csv'), 
            row.names = FALSE)
}

4.1 Visualize and assess

4.1.1 First draw of collapse free generated ensembles

tol archive draw n_stitched mean_gen var_gen time_period mean_target var_target bias E1 E2
0.1 with_ssp245 1 25 0.0650757 1.5313229 hist+future 0.0440780 1.5346484 0.0209977 0.0136824 0.9978331
0.1 without_ssp245 1 10 0.0381270 1.5216680 hist+future 0.0440780 1.5346484 0.0059510 0.0038778 0.9915418
0.1 with_ssp245 1 25 -0.9558477 0.4287856 hist -0.9783103 0.4366893 0.0224627 0.0514386 0.9819009
0.1 without_ssp245 1 10 -0.9761643 0.4356694 hist -0.9783103 0.4366893 0.0021460 0.0049143 0.9976645
0.1 with_ssp245 1 25 2.0238240 0.8087675 future 2.0056370 0.8086868 0.0181870 0.0224895 1.0000997
0.1 without_ssp245 1 10 1.9841510 0.7944827 future 2.0056370 0.8086868 0.0214860 0.0265689 0.9824355

In general, all of our emulations have E1 close to 0 and E2 close to 1. As always, the question is sort of ‘how close is close enough’.

Again, it’s good to remember tgat because the draw of generated ensemble Tgavs is random, there’s no guarantee that for an individual draw, using the archive including the SSP245 points will guarantee better fits. It will certainly create a larger potential space of better fits to draw from. For SSP585, I suspect this will be more of a factor.

4.1.2 And the second:

tol archive draw n_stitched mean_gen var_gen time_period mean_target var_target bias E1 E2
0.1 with_ssp245 10 24 0.0647906 1.5320185 hist+future 0.0440780 1.5346484 0.0207126 0.0134966 0.9982863
0.1 without_ssp245 10 11 0.0391560 1.5168975 hist+future 0.0440780 1.5346484 0.0049220 0.0032072 0.9884333
0.1 with_ssp245 10 24 -0.9569903 0.4285913 hist -0.9783103 0.4366893 0.0213200 0.0488220 0.9814558
0.1 without_ssp245 10 11 -0.9709116 0.4318697 hist -0.9783103 0.4366893 0.0073987 0.0169427 0.9889632
0.1 with_ssp245 10 24 2.0251842 0.8067235 future 2.0056370 0.8086868 0.0195472 0.0241715 0.9975722
0.1 without_ssp245 10 11 1.9770765 0.8020489 future 2.0056370 0.8086868 0.0285605 0.0353171 0.9917916

4.1.3 Aggregate metrics

We want to know the average errors across many draws so that we can get a sense of the error of our method and how that relates to hyperperameters:

# read in metrics for every draw
files <- list.files(path = paste0('generated_ensembles/', esm_name),
                    pattern = paste0('metrics_', esm_experiment), 
                    full.names = TRUE) 

generated_metrics <- list()
for (i in 1:length(files)){
  generated_metrics[[i]] <- read.csv(files[i], stringsAsFactors = FALSE)
}
generated_metrics <- do.call(rbind, generated_metrics)

generated_metrics %>% 
  group_by(tol, archive, time_period) %>%
  summarize(aggregate_E1 = mean(E1),
            aggregate_E2 = mean(E2)) %>%
  ungroup %>%
  mutate(total_draws = length(unique(generated_metrics$draw))) ->
  aggregate_metrics

kable(aggregate_metrics)
tol archive time_period aggregate_E1 aggregate_E2 total_draws
0.1 with_ssp245 future 0.0196261 0.9960188 10
0.1 with_ssp245 hist 0.0435020 0.9848697 10
0.1 with_ssp245 hist+future 0.0116808 0.9979305 10
0.1 without_ssp245 future 0.0399868 0.9963976 10
0.1 without_ssp245 hist 0.0147206 1.0023802 10
0.1 without_ssp245 hist+future 0.0085888 0.9915438 10
ggplot(data = aggregate_metrics, aes(x = aggregate_E1, y = aggregate_E2, color = factor(tol), shape = archive)) +
  geom_point() +
  facet_wrap(~time_period, nrow =3) +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = 0)

5 Validate - window jumps

Pull off the years that end each window and that start each window; these are going to be years we pay particular attention to.

window_start_years <- unique(archive_subset1$start_yr)
window_end_years <- unique(archive_subset1$end_yr)

We are particularly concerned about the ‘jumps’ at each of these years being ‘wrong’. And the idea is that they would be wrong relative to both the target data and, I think, the generated ensembles in between the jump years.

So let’s work with the jumps themselves for validating. Take the year to year diff within each realization (target or generated ensemble member). We’ll look at the diff and abs(diff). We’ll visualize the jumps with a scatter view and then get into more quantitative metrics. I spent a while going back and forth between starting from fundamentals like this and tracking down a package that looks at more sophisticated methods. As I started reading through packages, I ended up deciding to stick to fundamentals because they’ll translate from R to python better and will provide a more solid understanding of what’s going on before we look at anything more sophisticated.

5.1 calculate jumps

# read in all generated ensembles.
files <- list.files(path = paste0('generated_ensembles/', esm_name), 
                    pattern = paste0('generated_ensembles_', esm_experiment),  
                    full.names = TRUE) 

generated_tgavs <- list()
for (i in 1:length(files)){
  generated_tgavs[[i]] <- read.csv(files[i], stringsAsFactors = FALSE)
}
generated_tgavs <- do.call(rbind, generated_tgavs)


# make a helper function for calculating all the jumps because you have Nyears-1 jumps and
# dply functions don't like that. 
get_jumps <- function(data){
  data.frame(left_year = data$year[1:(length(data$year)-1)],
             jump = diff(data$value),
             abs_jump = abs(diff(data$value))) %>%
    mutate(variable = unique(data$variable),
           stitching_id = unique(data$stitching_id),
            tol = unique(data$tol),
           archive = unique(data$archive),
           draw = unique(data$draw),
           jump_year = if_else(left_year %in% window_end_years, 'window_transition', 'in_window')) ->
    x
return(x)
}

# calculate the yearly jumps
generated_list <- split(generated_tgavs, 
                        list(generated_tgavs$tol,
                        generated_tgavs$archive,
                        generated_tgavs$draw, 
                        generated_tgavs$stitching_id),
                        drop=TRUE)

lapply(generated_list, get_jumps) %>%
  do.call(rbind, .) ->
  generated_jumps

target_list <- split(original_data,
                     list(original_data$model, 
                          original_data$experiment,
                          original_data$ensemble),
                     drop = TRUE)

lapply(target_list, get_jumps) %>% 
  do.call(rbind, .) ->
  target_jumps

5.2 Visualize jumps

Generated jump year values during window transitions vs all target jump year values:

ggplot() + 
  geom_point(target_jumps, mapping = aes(x = left_year, y = jump), alpha = 0.5) +
  geom_point(filter(generated_jumps,
                    jump_year == 'window_transition'), mapping = aes(x = left_year, y = jump), alpha = 0.1, color='red') +
    facet_grid(tol~archive)

ggplot() + 
  geom_point(filter(target_jumps,
                    jump_year == 'window_transition'), mapping = aes(x = left_year, y = jump), alpha = 0.5) +
  geom_point(filter(generated_jumps,
                    jump_year == 'window_transition'), mapping = aes(x = left_year-1, y = jump), alpha = 0.1, color='red' ) +
    facet_grid(archive~tol) +
  ggtitle('note the years are staggered just for visual clarity')

6 distributions of jumps

What’s the appropriate comparison though?

  1. All years target (250 jumps X 25 ensemble members) vs all years generated (250 jumps X 5-50 ensemble members depending on tol+archive X Number draws under consideration)
  2. All years target (250 jumps X 25 ensemble members) vs window transition years generated (28 jumps X 5-50 ensemble members depending on tol+archive X Number draws under consideration)
  3. Window transition years target (28 jumps X 25 ensemble members)vs window transition years generated (28 jumps X 5-50 ensemble members depending on tol+archive X Number draws under consideration)

From the perspective of the target data, there’s nothing special about those transition years, which basically knocks out 3. I’ll plot it for completeness but I definitely don’t think it’s the right comparison

6.1 Full set of draws for method validation

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated'), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years, full ensemble (target or generated), \nTotal number of draws for generated ensembles ', length(unique(generated_jumps$draw))))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition'), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years target ensemble, \nWindow years generated ensemble, \nTotal number of draws for generated ensembles ', length(unique(generated_jumps$draw))))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target')%>% filter(jump_year == 'window_transition'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition'), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('Window years target ensemble, \nWindow years generated ensemble, \nTotal number of draws for generated ensembles ', length(unique(generated_jumps$draw))))

6.2 Single draw for use case validation

  • for the couple draws we look at, all the findings are qualitatively the same as the above yay.
focus_draw <- 1

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years, full ensemble (target or generated), \nFocusing on a single draw of generated ensemble ', focus_draw))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target')%>% filter(jump_year == 'window_transition'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('Window years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))

focus_draw <- 2

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years, full ensemble (target or generated), \nFocusing on a single draw of generated ensemble ', focus_draw))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target')%>% filter(jump_year == 'window_transition'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('Window years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))

6.3 Apply E1, E2 to jumps

Effectively quantifying how similar vs not the above distributions are, since they’re all pretty normal looking. I suppose we could do a KS test to support further if needed?

So we will compare all the years of generated data to all the years of target data, and the window transition years of generated data to all the years of target data. Individual draw and aggregate across draws errors

6.3.1 All the years for generated data

Individual draws:

# values of the target data
mean_target <- mean(target_jumps$jump)

var_target <- sd(target_jumps$jump)

  
for (drawID in 1:Ndraws){
  # Metrics for this draw #################################################################
  generated_jumps %>%
    filter(draw == drawID) ->
    generated_ensemble_draw 
  
  # Values of the generated data and then the metrics
  generated_ensemble_draw %>% 
    group_by(tol, archive, draw) %>%
    summarize(mean_gen = mean(jump), 
              var_gen = sd(jump)) %>%
    ungroup %>%
    mutate(mean_target = mean_target,
           var_target = var_target,
           bias = abs(mean_target - mean_gen),
           E1 = bias/var_target,
           E2 = var_gen/var_target) ->
    metrics
  
  
  
  generated_ensemble_draw %>%
    select(stitching_id, tol, archive, draw) %>%
    distinct %>%
    group_by(tol, archive, draw) %>%
    summarise(n_stitched = n()) %>%
    ungroup  %>%
    left_join(metrics, by = c('tol', 'archive', 'draw')) ->
    plotting_tbl
  
   write.csv(plotting_tbl %>%
              mutate(draw = drawID),
            paste0('generated_ensembles/', esm_name ,'/jump_errors_allyears_', esm_experiment, '_draw', drawID, '.csv'), 
            row.names = FALSE)
}


# Plot the first draw's errors:
plotting_tbl <- read.csv(paste0('generated_ensembles/',esm_name , '/jump_errors_allyears_',esm_experiment,'_draw1.csv'),  
                         stringsAsFactors = FALSE)

kable(plotting_tbl)
tol archive draw n_stitched mean_gen var_gen mean_target var_target bias E1 E2
0.1 with_ssp245 1 25 0.0176517 0.1216254 0.0177774 0.1176741 0.0001257 0.0010682 1.033578
0.1 without_ssp245 1 10 0.0180020 0.1246752 0.0177774 0.1176741 0.0002246 0.0019085 1.059495
ggplot(data = plotting_tbl, aes(x = E1, y = E2, color = factor(tol), shape = archive)) +
  geom_point() +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = 0) + 
  ggtitle('Error metrics for single Draw 1, comparing all \nyears of generated jumps to all years of target jumps')

Aggregated across draws:

# read in metrics for every draw
files <- list.files(path = paste0('generated_ensembles/', esm_name),
                    pattern = paste0('jump_errors_allyears_', esm_experiment), 
                    full.names = TRUE) 

generated_metrics <- list()
for (i in 1:length(files)){
  generated_metrics[[i]] <- read.csv(files[i], stringsAsFactors = FALSE)
}
generated_metrics <- do.call(rbind, generated_metrics)

generated_metrics %>% 
  group_by(tol, archive) %>%
  summarize(aggregate_E1 = mean(E1),
            aggregate_E2 = mean(E2)) %>%
  ungroup %>%
  mutate(total_draws = length(unique(generated_metrics$draw))) ->
  aggregate_metrics

kable(aggregate_metrics)
tol archive aggregate_E1 aggregate_E2 total_draws
0.1 with_ssp245 0.0021871 1.042077 10
0.1 without_ssp245 0.0014537 1.066417 10
ggplot(data = aggregate_metrics, aes(x = aggregate_E1, y = aggregate_E2, color = factor(tol), shape = archive)) +
  geom_point() +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = 0) + 
  ggtitle(paste('Aggregate errors (jumps for all years), aggregated across Ndraws ', unique(aggregate_metrics$total_draws) ))

6.3.2 Window transition years for generated data

Individual draw:

# values of the target data
mean_target <- mean(target_jumps$jump)

var_target <- sd(target_jumps$jump)

  
for (drawID in 1:Ndraws){
  # Metrics for this draw #################################################################
  generated_jumps %>%
    filter(draw == drawID,
           jump_year=='window_transition') ->
    generated_ensemble_draw 
  
  # Values of the generated data and then the metrics
  generated_ensemble_draw %>% 
    group_by(tol, archive, draw) %>%
    summarize(mean_gen = mean(jump), 
              var_gen = sd(jump)) %>%
    ungroup %>%
    mutate(mean_target = mean_target,
           var_target = var_target,
           bias = abs(mean_target - mean_gen),
           E1 = bias/var_target,
           E2 = var_gen/var_target) ->
    metrics
  
  
  
  generated_ensemble_draw %>%
    select(stitching_id, tol, archive, draw) %>%
    distinct %>%
    group_by(tol, archive, draw) %>%
    summarise(n_stitched = n()) %>%
    ungroup  %>%
    left_join(metrics, by = c('tol', 'archive', 'draw')) ->
    plotting_tbl
  
   write.csv(plotting_tbl %>%
              mutate(draw = drawID),
            paste0('generated_ensembles/', esm_name ,'/jump_errors_transitionyears_',esm_experiment,'_draw', drawID, '.csv'), 
            row.names = FALSE)
}


# Plot the first draw's errors:
plotting_tbl <- read.csv(paste0('generated_ensembles/', esm_name ,'/jump_errors_transitionyears_',esm_experiment,'_draw1.csv'),  
                         stringsAsFactors = FALSE)

kable(plotting_tbl)
tol archive draw n_stitched mean_gen var_gen mean_target var_target bias E1 E2
0.1 with_ssp245 1 25 0.0147275 0.1576160 0.0177774 0.1176741 0.0030499 0.0259185 1.339428
0.1 without_ssp245 1 10 -0.0002835 0.1681503 0.0177774 0.1176741 0.0180609 0.1534825 1.428949
ggplot(data = plotting_tbl, aes(x = E1, y = E2, color = factor(tol), shape = archive)) +
  geom_point() +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = 0) + 
  ggtitle('Error metrics for single Draw 1, comparing window transition years of \ngenerated jumps to all years of target jumps')

Aggregated across draws:

# read in metrics for every draw
files <- list.files(path = paste0('generated_ensembles/', esm_name), 
                    pattern = paste0('jump_errors_transitionyears_',esm_experiment ), 
                    full.names = TRUE) 

generated_metrics <- list()
for (i in 1:length(files)){
  generated_metrics[[i]] <- read.csv(files[i], stringsAsFactors = FALSE)
}
generated_metrics <- do.call(rbind, generated_metrics)

generated_metrics %>% 
  group_by(tol, archive) %>%
  summarize(aggregate_E1 = mean(E1),
            aggregate_E2 = mean(E2)) %>%
  ungroup %>%
  mutate(total_draws = length(unique(generated_metrics$draw))) ->
  aggregate_metrics

kable(aggregate_metrics)
tol archive aggregate_E1 aggregate_E2 total_draws
0.1 with_ssp245 0.0323019 1.339463 10
0.1 without_ssp245 0.0565644 1.470436 10
ggplot(data = aggregate_metrics, aes(x = aggregate_E1, y = aggregate_E2, color = factor(tol), shape = archive)) +
  geom_point() +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = 0) + 
  ggtitle(paste('Aggregate errors comparing window transition years of \ngenerated jumps to all years of target jumps \naggregated across Ndraws ', unique(aggregate_metrics$total_draws) ))